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SECTION  I 
INTRODUCTION 


During  the  past  several  years  the  authors  have  generated  analytical 
techniques  for  predicting  the  mutual  aerodynamic  Interference  between 
externally  carried  stores  and  between  the  stores  and  the  aircraft  (Ref- 
erences 1 through  10) . Although  some  solution  techniques  using  surface 
singularities  are  considered,  the  basic  approach  utilizes  the  Ranklne 
method  of  superposition  of  an  axial  source-sink  distribution  to  form  a 
closed  body.  The  interference  flow  fields  are  then  superimposed  on  the 
Isolated  solutions  for  the  external  stores  (finned  bodies  of  revolution) 
and  the  aircraft  components  (fuselage,  wing,  and  pylon). 

The  work  previously  performed  In  this  general  area  has  been  limited 
to  solutions  of  the  Prandtl-Glauert  equation,  a form  of  the  governing 
differential  equation  which  was  linearized  by  assuming  Invlscld  fluid, 
small  velocity  perturbations,  and  avoiding  the  transonic  regime.  Com- 
parison of  this  analytical  work  with  wind  tunnel  test  data  has  shown, 
however,  that  the  linearized  solutions  will  not  predict  accurately  the 
compressible  flow  effects  at  conditions  approaching  the  critical  mach 
number.  Further,  In  the  case  of  the  triple  ejector  rack  (TER)  with  stores, 
the  pressures  observed  experimentally  are  indicative  of  large  velocity 
perturbations,  l.e.,  near  stagnation  values  forward  of  the  bomb  rack.  It 
Is  felt,  therefore,  that  both  transonic  and  large  perturbations  terms  will 
have  to  be  considered. 


In  addition  to  the  large  velocity  perturbation  and  transonic  terms. 

It  Is  very  likely  that  the  viscous  effects  will  have  to  be  Included  In 
order  to  predict  the  loading  for  some  configurations.  This  is  particular- 
ly true  In  the  case  of  multiple  weapons  carriage  such  as  with  TERs  and 
multiple  ejector  racks  (MERs).  For  example.  In  the  case  of  three  M117 
bombs  on  a TER,  the  boundary  layers  between  bombs  coalesce.  Thus,  the 
usual  simple  boundary  layer  displacement  corrections  cannot  be  used.  Also, 
In  the  case  of  a shockwave-boundary  layer  interaction,  the  solutions  for 
the  boundary  layer  and  the  potential  flow  are  mutually  dependent;  there- 
fore, the  significant  terms  in  the  governing  differential  equations  are 
nonlinear. 


While  it  Is  quite  desirable  to  obtain  a solution  which  Includes  all 
of  the  nonlinear  terms,  more  restrictive  approaches  are  generally  more 
fruitful.  Because  of  this,  several  approximate  as  well  as  exact  approaches 
were  simultaneously  pursued  In  this  effort. 


At  the  time  that  this  effort  was  Initiated,  it  was  felt  that  the 
approach  most  likely  to  succeed  for  the  Invlscld  flow  case  would  be  to 
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generate  higher  order  effects  based  on  the  current  Image-systen  solution. 
Basically,  the  approach  would  be  to  evaluate  the  effects  of  the  nonlinear 
terns  which  previously  have  been  neglected  from  the  governing  differential 
equation  for  Invlscld  flow.  This  would  be  done  by  the  successive  approx- 
imation technique  or  by  the  more  orderly  method  of  asymptotic  expansion. 
Considerable  effort  was  put  Into  this  approach,  and  the  results  are 
discussed  In  this  report. 

Initially  it  was  felt  that  the  best  approach  for  the  viscous  flow 
problem  would  be  to  use  the  current  potential  flow  solution  (see  References 
1 through  10)  for  Che  far-fleld  boundary  conditions  and  use  an  Imbedded 
numerical  solution  technique  for  the  boundary  layer  and  boundary  layer- 
shock  wave  solutions.  This  concept  was  pursued  using  classical  nume  ical 
marching  techniques.  It  became  evident,  however,  that  a more  general 
numerical  solution  Is  required. 

In  addition,  other  numerical  approaches  have  been  studied  which,  along 
with  the  efforts  mentioned  above,  are  discussed  In  detail  In  the  following 
sections  of  this  report. 
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SECTION  II 


APPROXIMATE  SOLUTIONS 


INTRODUCTION 

The  current  solutions  given  in  References  1 through  10  are  limited 
to  the  case  of  slender  bodies  for  which  the  flow  is  all  either  subsonic 
or  supersonic.  Specifically,  this  previous  work  represents  solutions  to 
the  Pcandtl-Glauert  equation  which  is  the  linearized  velocity  potential 
equation  given  as  Equation  (1) . 


(l-M-^)((>  + 

“>  XX 


yy 


zz 


(1) 


Noting  that  Equation  (1)  becomes  the  Laplace  equation  if  the  free 
stream  mach  number  is  zero,  it  is  seen  that  for  this  case  the  restriction 
of  small  velocity  perturbations  can  be  removed.  This  is  the  so-called 
zeroth  solution  and  can  be  used  as  a basis  for  generating  the  compres- 
sibility effects  at  higher  mach  numbers  by  a successive  approximation 
technique,  or  by  the  more  orderly  method  of  assuming  a solution  in  the 
form  of  a series.  A classical  series  solution  known  as  the  Rayleigh- 
Janzen  method  has  been  used  in  a nearly  exact  manner  for  a few  simple 
problems  as  noted  by  Shapiro  (Reference  11) . Although  this  type  of  solu- 
tion technique  has  been  greatly  extended  (see  Reference  12)  by  other 
investigators,  it  has  been  applied  only  to  simple  isolated  bodies.  It 
was  felt,  however,  that  because  of  the  current,  closed  form,  analytical 
solutions  for  aerodynamic  Interference  (see  References  1 through  10)  the 
Rayleigh-Janzen  method  could  be  used  for  this  problem.  In  order  to 
verify  this  approach,  several  formulations  of  the  solution  were  generated 
for  an  Isolated  slender  body  as  follows: 

SLENDER  BODY  SOLUTION 

The  usual  form  of  the  velocity  potential  equation  is  written  as 
follows  (for  two  dimensions): 

(ip2  - a^)  Ip  + 2ip  p p + (p^  - a^)  p - 0 (2) 

X XX  X y xy  y yy 

where  a is  the  sonic  speed  and  can  be  expressed  as 

a2  - a2  [1  + (^)  m2  (v2  - p2  _ (3) 
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Equation  (2)  Is  the  form  of  the  velocity  potential  equation  usually 
found  In  the  literature  as  the  most  general  form  (for  tvro-dlmenslons) . 
For  this  study,  however.  It  was  found  that  a more  useful  form  Is  given 
by  Equation  (4)  or  its  equivalent  Equation  (5) . 


7(ji*Va2 

724,  . - (-^)  7(|)*71n  a2 


which  hold  for  both  two-  and  three-dimensional  cases. 

Expanding  Equation  (4)  in  terms  of  the  perturbation  velocity 
potential  yields  a nondlmenslonal  form  (3  dimensions) : 

[1  - (^)  m2  (2^^  + + ^2  + ^2)]  ^2^  , 

m2 

-T—  [{41  +1)  T — + T — + ♦ X — 1 (2<^  + 4>^  •l*^)  (6) 

2 3x  y 3y  z 3y  x x y z 

Since  we  are  looking  for  the  compressibility  effects,  we  can  set 

e ■ m2,  where  e will  be  our  perturbation  parameter,  and  assume  a solution 
of  the  form 

(Kx,y,z;e)  - '('^^^x.y.z)  + e 4>  ^^\x,y,z)  + e2  4i^^^(x,y,z)  + ...  (7) 

The  assumed  solution  can  be  substituted  In  Equation  (6)  and  (in 
principle)  each  successive  order  of  magnitude  term  can  be  solved.  For 
small  perturbations,  however.  Equation  (6)  can  be  simplified  by  neglecting 
velocity-squared  terms  in  comparison  to  first-order  terms  and  after  some 
manipulation  reduces  to 


[1  - (y-D  m2  4>  J a2724,  , j,2^ 

“>  X * XX 


For  slender  axisymmetric  bodies  at  zero  angle  of  attack  the  body 
boundary  conditions  can  be  expressed  as 


♦y  cos  0 + 4I2  ® 
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Substituting  Equation  (7)  into  Equations  (8)  and  (9),  and  equating  terms  In 
powers  of  e gives  a set  of  order  of  magnitude  equations.  The  first  three, 
zeroth,  first,  and  second  orders,  are  given  In  Appendix  A. 


The  numerical  results  of  applying  this  solution  to  an  M117  bomb  shape 
are  presented  as  Figures  1 and  2.  On  Figure  1,  the  theoretical  results 
through  the  second-order  terms  are  compared  with  the  wind  tunnel  test  results 
of  Reference  13.  Figure  2 shows  a comparison  with  experiment  for  the  second- 

order  solution;  there  is  very  good  agreement  for  mach  numbers  0.5  and  0.8. 

At  mach  number  0.8,  however,  the  critical  free  stream  mach  number  has  been 
exceeded  as  evidenced  by  a slight  jump  In  the  pressure  distribution  at  an 
axial  station  of  about  X/D  ■ 1.5.  A much  larger  pressure  Jump  can  be  seen 

for  the  mach  number  case  of  M = 0.9.  The  pressure  on  the  nose  of  the  body  is 

predicted  quite  well,  however,  up  to  a value  of  X/D  - 1.0  where  the 
local  flow  has  become  supersonic. 

In  principle,  this  solution  can  be  applied  to  the  multibody  interference 
solution  as  outlined  In  References  3,  6,  8,  and  9.  In  the  case  of  three 
large  bombs  mounted  on  a triple  ejector  rack,  however,  there  will  be  nearly 
stagnation  pressures  in  the  vicinity  of  the  bomb  racks,  and,  therefore,  the 
small  perturbation  assumption  will  be  violated. 

It  was  felt  that  a more  exact  solution  would  be  needed  In  order  to 
obtain  a valid  solution  In  the  vicinity  of  stagnation  points;  therefore,  a 
modified  slender  body  solution  has  been  formulated. 

MODIFIED  SLENDER  BODY  SOLUTION 

A convenient  form  of  the  differential  equation  which  Is  valid  at  stag- 
nation points  can  be  obtained  by  using  Equatlon(5)  which,  written  In  terms 
of  the  non-dimensional  perturbation  velocity  potential,  becomes 

^ + V(fr-7)an  a2  (10) 

where  a^  Is  the  square  of  the  nondlmenslonal  speed  of  sound  and  In  terms  of 
the  perturbation  potential  function  can  be  expressed  as 

a2  - 1 _ m2(2({i  + + ^2  + ^2)  (12) 

2 “ X X y z 


For  adiabatic  flow,  the  non-dimensional  speed  of  sound  can  be  expressed  in 
terms  of  the  local  and  free  stream  mach  numbers  and  Is  given  by 


From  this  expression  It  can  be  seen  that  the  minimum  value  of  a^  for  sub- 
crltlcal  speeds  Is  given  for  M"0  and  This  gives  a value  of  a2*0.8333. 

Thus,  an  appropriate  expansion  for  In  Is  given  by  the  series. 


In  a2  . (a^-l)  - i(a2-l)2  + j(a2-l)3  +. . . (13) 

The  maximum  error  In  approximating  In  a^  with  the  first  term  of  the 
series  Is  8.5  percent,  and  if  the  first  two  terms  are  used,  the  error  is 
only  0.82  percent. 

It  follows,  therefore,  that  to  good  approximation.  Equation  (10)  can 
be  expressed  as 

+ V<>.V][(a2-l)  - i(a2-l)2]  (14) 

Substituting  Equation  (11)  into  Equation  (14)  yields 

l|^  + V^.VH(2^^+*2-H,2+^2)  + 

(■^) Y-(2i>^+<>2+()|2+,J,2 ) 2 ] (2^5 ) 

Examining  the  differential  operator  on  the  right-hand  side.  It  can  be 
seen  that  the  coefficients  vanish  at  stagnation  points  as  described  below: 

The  differential  operator  In  expanded  form  Is 


(l+(j)  ) — — + "7 — + ^ -7 — 

X 3x  y 3y  3z 


(16) 


At  stagnation  points  the  perturbation  velocities  are  given  by 


♦ - 0 (17) 

s 

Thus,  the  right-hand  side  of  Equation  (15)  Is  zero,  and  the  governing 
equation  at  stagnation  points  Is 
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72((.  - 0 


(18) 
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NOTE:  Equation  (18)  holds  at  stagnation  points  whether  the  term, 
tn  a^.  Is  or  Is  not  approximated  with  the  truncated  series. 

It  follows  that  the  term  on  the  right-hand  side  of  Equation  (15)  has 
negligible  effect  In  the  vicinity  of  stagnation  points. 

Now,  If  It  Is  assumed  that  the  velocity  perturbations  are  small  at 
points  a'.!.'<^y  from  stagnation  conditions.  It  can  be  assumed  that 

« 2(>  (19) 

X y z X 

To  good  approximation,  therefore,  Equation  (15)  can  be  expressed  as 
« l|^  + V.^.V]Im2  + (1^)  ^2)  (20) 


Equation  (20)  properly  accounts  for  stagnation  conditions,  small  velocity 
perturbations  away  from  stagnation  points,  and  conditions  up  to  the  critical 
mach  number.  The  Rayleigh-Janzen  solution  for  this  equation  Is  developed 
for  the  first-order  terms  In  Appendix  B.  Comparisons  of  this  result  with 
experiment  and  other  analytical  results  are  shown  In  Figures  3 and  4. 

Before  proceeding  to  a more  complicated  geometry  or  analysis,  however. 

It  Is  of  Interest  to  note  the  following: 

1.  Comparison  of  the  Rayleigh-Janzen  results  (this  analysis,  both 
slender  body  and  modified  slender  body)  with  the  solution  using  Gothert's 
Rule  Indicates  that  little  Improvement  can  be  expected  by  going  to  this  more 
complicated  technique.  This  certainly  seems  true  at  least  for  bodies  of 
revolution  at  zero  angles  of  attack  and  subcrltlcal  mach  numbers  as  Indicated 
by  Figures  3 and  4. 

2.  From  the  discussion  leading  to  Equation  (18),  the  Laplace  equation 
holds  in  the  near  vicinity  of  stagnation  points.  This  Implies  that  the 
small  perturbation  solutions  are  not  significantly  In  error  near  stagnation 
points.  This  conclusion  Is  borne  out  by  the  agreement  between  the  theories 
and  experiments  shown  on  Figures  1 through  4. 

3.  The  analytical  approaches  discussed  above  are  valid  only  for  pure 
subsonic  flow.  This  can  be  aeen  on  Figure  2.  There  is  only  a alight  jump 
In  the  experimental  pressure  coefficient  at  M^*0.8,  but  at  H ■0.9,  thla 
Jump  la  quite  prominent.  In  fact,  the  dlacrepancy  between  the  analytical 
and  experimental  pressure  coefficients  la  quite  large  for  the  M^"0.9  case. 
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- RAYLEKJH-JANZEN  (FIRST  ORDER,  MODIFIED  SMALL 
DISTURBANCE) 

RAYLEIGH-JANZEN  (SECOND  ORDER,  SMALL 

DISTURBANCE) 

— GOTHERT'S  RULE  (SMALL  DISTURBANCE) 
EXPERIMENT  (REFERENCE  13) 


oO^ 


5 X/D  6 


(A).  Pressure  Distribution  M * 0.8 


MODEL  SHAPE 


(B).  Schematic  of  Body  Profile 


Figure  3.  Coaparieon  of  Raylelgh-Janzen  Solution  and  Solution  Using 
Gothert's  Rule,  M117  Bomb  Shape,  a ■ 0“,  M ■ 0.8 
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It  seems  evident,  therefore,  that  the  Rayleigh-Janzen  solution  tech- 
^ nique,  at  least  through  the  second-order  solution,  will  not  be  significantly 

better  at  predicting  the  aerodynamic  loads  near  and  at  the  critical  mach 
niimKor  than  the  much  simpler  technique  using  Gothert's  Rule.  It  still  may 
be  possible,  however,  to  use  the  Rayleigh-Janzen  solution  at  the  critical 
mach  nu]id>er  if  the  higher  order  terms  of  the  differential  equation  can  be 
retained  and  if  the  higher  order  solutions  can  be  generated.  An  attempt  to 
retain  the  higher  order  effects  is  proposed  as  the  next  approach  to  the 
problem. 

SUCCESSIVE  NUMERICAL  APPROXIMATION 

The  actual  analytical  development  of  the  higher  order  solutions  by  the  • 

Rayleigh-Janzen  method  is  a horrendous  Job  even  for  the  most  simple  case. 

Even  for  a single  isolated  body  for  which  the  zeroth  solution  is  obtained  by 
determining  the  centerline  source  distribution,  the  third-order  solution  is 
most  unwieldy.  It  is  proposed,  therefore,  that  a numerical  solution  be 
attempted  as  follows: 

First,  it  is  assumed  that  the  zeroth  solution  has  been  obtained.  This 
solution  is  not  limited  to  centerline  sources  and  doublets,  but  can  be 
applied  to  any  form  of  solution,  including  source  panels,  vortex  panels,  and 
vortex  sheets.  Assuming  the  zeroth  solutions  and  its  derivatives  are  known, 
a higher  order  solution  can  be  generated  as  outlined  in  the  following  steps. 

1.  Write  Equation  (6)  as 

* 2 

m2 

V2^  - + V(|)*V][24.  + (21) 

2a2  3*  ^ 

where  I 

a2  - 1 - (^)  m2  [2^^  + (22)  i 

and  ^ is  the  non-dimensional  velocity  potential  for  uniform  flow  in  the 
x-directlon. 

If  the  usual  Rayleigh-Janzen  method  were  followed.  Equation  (7)  would  j 

be  substituted  into  Equations  (21)  and  (22)  and  the  coefficients  of  powers 
of  the  mach  nunfcer  equated.  This,  of  course,  would  result  in  the  ordered 
solutions  as  previously  considered.  It  is  easily  seen  that  the  zeroth 
solution  would  be  the  solution  obtained  by  solving  the  Laplace  equation 
written  as 


The  solution  of  this  equation  will  also  be  taken  as  the  first  step  In  this 
method.  However,  the  actual  method  of  solving  Equation  (23)  Is  unimportant 
to  this  approach. 

2.  The  next  (and  all  successive)  steps  will  follow  the  usual  successive 
approximation  method.  Thus,  to  the  next  approximation.  Equation  (21)  can  be 
written  as 


where  the  complete  boundary  condition  Is  to  be  applied  to  the  solution  for 
and  F(^^^^)  Is  the  right-hand  side  of  Equation  (21)  expressed  In 
terms  of  . 

3.  Write  the  first  approximation  as  the  sum  of  the  particular  and 
complimentary  solutions  which  can  be  expressed  as 

♦ “ + ♦p  (25) 

where  the  particular  solution  must  satisfy  the  partial  differential  equation 

(26) 

while  the  complimentary  solution  must  satisfy  the  equation 

- 0 (27) 


and  Equation  (25)  must  satisfy  the  boundary  conditions. 

The  particular  solution  can  be  determined  throughout  the  flow  field  by 
numerical  techniques.  However,  If  the  only  points  of  Interest  for  a single 
body  are  those  lying  on  the  body  surface,  the  size  of  the  numerical  grid  can 
be  greatly  reduced.  If  two  or  more  bodies  with  mutually  Interfering  flow 
fields  are  considered,  a grid  system  which  envelops  all  bodies  must  be  used. 
This  system  could  still  be  two  or  more  orders  of  magnitude  smaller  than  that 
required  by  other  numerical  techniques. 

4.  Successive  approximations  can  be  generated  by  replacing  the  zeroth 

solution,  with  the  first  approximation  and  repeating  steps  2 

and  3. 
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The  classical  Raylelgh-Janzen  method  is  valid  only  for  subcrltlcal  mach 
numbers.  Conceptually,  however,  this  approach  could  be  made  useful  for 
transonic  and  perhaps  supersonic  flows  by  taking  Into  account  von  Kirmdn's 
Rule  for  Forbidden  Signals. 

An  alternate  formulation  for  this  method  which  might  greatly  speed  up 
convergence  for  the  subsonic  case  and  help  to  Implement  the  supersonic 
solution  Is  as  follows: 

1.  Write  Equation  (21)  as 

,2,  . [|_*  „.vH2»__+,2*^2+,2,  . ,28) 

2.  For  the  subsonic  case,  this  results  In  the  zeroth  solution  which 
Is  the  solution  to  the  Prandtl-Glauert  equation 

(l-M^)  <(.  +4.  + ((>  - 0 (29) 

“>  XX  yy  zz 

The  solution  of  this  equation  using  Gothert's  Rule  Is  compared  with  the 
Raylelgh-Janzen  method  In  Figures  3 and  4. 

3.  For  the  supersonic  case  the  zeroth  solution  Is  the  solution  of  the 
classical  wave  equation  written  as 


(m2-1)  4> 

<*>  XX 


zz 


0 


(30) 


Several  techniques  are  known  for  solving  this  equation. 

Further  evaluation  of  the  use  of  Equations  (28),  (29),  and  (30)  will 
be  required  before  their  use  In  this  successive  approximation  scheme  can  be 
validated. 
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EXACT  NUMERICAL  SOLUTION,  PURVIS  METHOD 


INTRODUCTION 

The  analytical  approaches  found  in  References  1 through  10  and  the 
methods  discussed  In  the  previous  section  do  not  seem  to  be  adequate  at 
speeds  approaching  the  critical  mach  number  which  Is  the  speed  range  of 
most  Interest  for  this  study.  This  Inadequacy  seems  to  be  due  to  the  fact 
that  the  previous  methods  have  used  linearized  potential  flow  theory,  and 
at  speeds  approaching  critical  mach  number,  the  neglected  nonlinear  terms 
become  dominant. 

In  this  study,  an  attempt  Is  made  to  provide  an  analytical  foundation  for 
a method  which  may  ultimately  be  used  to  solve  the  full  nonlinear  potential 
flow  field  about  any  configuration.  The  problem  Is  doubly  complicated  due  to 
the  nonlinear  partial  differential  equation  (PDE)  Involved  and  the  fact  that 
the  domain  Is  fully  elliptic;  l.e.,  disturbances  at  any  point  are  felt  every- 
where In  the  flow  field.  Since  the  problem  Is  nonlinear,  it  Is  felt  that  a 
numerical  scheme  must  be  used. 

The  major  numerical  methods  for  solving  a nonlinear  PDE  of  this  type 
are  Finite  Difference  (Relaxation  and  Time  Dependent  asymptotic  approach  to 
steady  state).  Finite  Elements,  Method  of  Integral  Relations,  and  Method 
of  Projections. 

Comprehensive  surveys  and  descriptions  of  these  techniques  and  applica- 
tions may  be  found  In  references  such  as  References  14  and  15  and  will  not 
be  discussed  here.  Most  of  these  methods  have  been  eliminated  from  consider- 
ation due  to  excessive  computer  run  times  or  storage  requirements  and  coie- 
plexlty  of  the  coding. 

The  current  work  presents  a much  simpler  numerical  method  for  solving 
the  full  nonlinear  potential  equation  than  those  found  In  the  literature. 

The  ensuing  analysis  uses  some  of  the  more  desirable  concepts  of  several  of 
the  above  techniques  to  minimize  solution  time  and  complexity  and  simplify 
the  treatment  of  body  boundary  conditions 

In  general,  the  approach  blends  the  finite  element  and  finite  difference 
concepts  into  a method  for  solving  the  full  potential  equations  in  a conserva- 
tion form.  The  problem  formulation  leads  to  a large  system  of  equations 
which  is  solved  iteratively.  The  solution  does  not  require  storing  or  invert- 
ing a large  coefficient  matrix  and  proceeds  in  such  a manner  that  disturbances 
can  propagate  and  be  felt  at  all  points  in  the  field  during  the  solution. 
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GOVERNING  EQUATION 


Under  the  assumptions  of  steady,  adiabatic,  Inviscld  flow,  the  equations 
for  conservation  of  mass,  momentum,  and  energy  and  the  ideal  gas  equation  of 
state  may  be  combined  into  a single  differential  equation.  A further  re- 
striction of  irrotationality  allows  the  introduction  of  the  velocity  poten- 
tial. The  governing  equation  may  then  be  expressed  as 


V • (pV()))  ■ 0 


(31) 


where  is  the  nondimenslonalized  total  velocity  potential.  Since  the  flow 
is  Isentroplc,  the  non-dimensional  density  p may  be  expressed  in  terms  of 
the  velocity  potential  as 

1 

p - [1  + ^ m2  (1  - (32) 

Equation  (31) , when  expanded  using  the  components  of  the  del  operator  V and 
put  in  terms  of  the  non-dimensional  perturbation  potential,  is  the  familiar 
full  potential  equation. 

DIVERGENCE  INTEGRAL 

To  formulate  the  problem  for  a computer  solution,  the  divergence  theorem 
is  applied  to  the  Integral  form  of  the  governing  Equation  (31).  Integrating 
over  an  arbitrary  volume  V yields 


/ V*(pV(j>)dV  - ^ pV4.*ndS  « 0 (33) 

V S 

where  S is  the  boundary  surface  enclosing  V and  n is  the  unit  outward  normal 
on  S. 


The  surface  Integral  form  of  Equation  (33)  is  solved  numerically  as 
discussed  in  the  following  subsections. 

GENERAL  SOLUTION 

As  in  other  numerical  methods,  the  domain  in  which  Equation  (33)  is  to 
be  solved  is  subdivided  with  a rectangular  mesh  or  grid,  and  initial  values 
of  p arc  assumed  at  the  center  of  each  cell  in  the  grid.  The  surface 
Integral  form  of  Equation  (33)  is  applied  to  each  cell  in  the  grid  and  is 
written  in  terms  of  the  potential  at  the  center  of  each  cell.  The  resulting 
system  of  equations  is  solved  by  successive  over-relaxation  and  subject  to 
the  usual  boundary  conditions  of  no  flow  normal  to  body  surfaces  and  uniform 
flow  at  infinity.  The  values  of  ^ thus  found  are  then  differenced,  and  new 
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values  for  p In  each  cell  are  obtained  using  Equation  (32) . This  process  Is 
repeated  until  values  of  p In  each  cell  have  converged  to  within  a desired 
percentage  for  several  successive  cycles. 

CELL  FORHULATION 

the  surface  Integral  form  of 
Figure  5 , Is : 


- 0 (34)  ; 

To  obtain  Equation  (34) , the  sign  on  Equation  (33)  was  changed  and  the 
Integral  was  evaluated  assuming  a linear  distribution  of  pu  or  pv  along  each 
surface  S.*  The  values  In  Equation  (34)  are  then  the  values  at  the  center  of  ' 

each  surface,  or,  equivalently,  the  average  of  pu  or  pv  along  S.  I 

Expressed  In  terms  of  values  at  the  center  of  the  cell,  the  terms  In 
Equation  (34)  are 


For  two-dimensional  and  axlsymmetrlc  flows. 
Equation  (33)  as  applied  to  cell  (l,j),  shown  In 


^iVi  ■ P3V3  ■ ‘’4V4 


and 


2 ^^i,j  ^i-l,j^ 

(35) 

Vx 

(36) 

Equation  (36)  Is  the  standard  centered  difference  formula  and  Is  accurate 
to  order  (Ax)^.  The  areas  In  Equation  (34)  differ  for  2-D  and  axlsymmetrlc 
flow.  Referring  to  Figure  S,  the  areas  for  2-D  flow  are 


s,  - y-  - y-  (37) 

j+i  J 

S,  ■ X - X (38) 

^ ®1+1  *1 

where  a unit  depth  has  been  assumed.  For  the  3-D  axlsymmetrlc  case,  the 
areas  are 

Sj^  - iry£  - iry^ 


J+1  i 

^ *l-*-l  *1 


(39) 

(40) 


*Equlvalent  to  assuming  a second-order  distribution  of  ^ along  each  cell 
surface  in  the  Finite  Element  method. 
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Using  the  appropriate  forms  of  Equations  (35)  through  (40),  Equation 
(34)  can  be  rearranged  into 


■ » 


(41) 


where,  for  example. 


A1  - - (p^^j  + p^_^^j)  S^/2Ax 


(42) 


and 


A5  - - (A1  + A2  + A3  + A4) 


(43) 


In  matrix  formulation,  this  is  a banded,  diagonally  dominant  matrix,  and 
thus  an  Iterative  method  such  as  accelerated  Gauss-Eledel  (successive  over- 
relaxation)  will  always  converge.  Note  that  Equation  (34)  Is  exact  In  the 
limit  of  vanishing  cell  size  and  that  only  Ist  derivatives  of  the  dependent 
variable  are  required.  The  finite  difference  equations  used  for  u and  v are 
accurate  to  order  (Ax)^  and  (Ay)^,  respectively. 

BODY  BOUNDARY  CONDITIONS 

One  of  the  most  difficult  aspects  of  finite  difference  methods  to  date 
has  been  the  treatment  of  body  conditions,  while  the  opposite  Is  true  In  the 
finite-element  method  since  the  body  surface  forms  the  sides  of  some  elements. 
The  present  formulation  combines  the  simplicity  of  calculations  of  the  finite- 
difference  method  with  the  relatively  sl]iq>le  but  exact  finite-element  treat- 
ment of  body  boundary  conditions.  Consider  the  cell  cut  by  the  body  surface 
as  shown  In  Figure  6.  The  surface  Integral  for  this  cell  Is  around  the 
boundary  shown  by  the  double  sol^d  line.  Since  the  body  boundary  condition  of 
no  flow  normal  to  the  surface  (V*n«0)  requires  that  there  be  no  flow  through 
surface  AB,  then  the  equivalent  expression  of  Equation  (34)  Is 


Pj^Uj^Sj^  - P3UJS2  - - 0 (44) 

where  S3  Is  now  the  area  of  segment  BC.  Note  that  this  equation  Is  exact 
for  any  size  cell,  and  that  no  Interpolation,  contour  mapping,  or  other 
scheme  Is  required  as  with  general  finite-difference  methods. 

STRETCHED  COORDINATES 

Since  It  Is  desirable  to  keep  the  number  of  cells  at  a minimum,  for  some 
problems  It  is  necessary  to  employ  nonlinear  cell  stretching  to  place  the 


19 


r 


outer  boundaries  at  infinity.  A stretched  coordinate  capability  is  provided 
using 


tan  C 


(45) 


tan  6 


(46) 


From  the  chain  rule,  the  finite-difference  derivatives  become,  for  example: 


= it  lii 

3x 


-1.1  J 

AC 


1 + 


(47) 


Hence,  a uniform  rectangular  grid  in  C,6  can  be  used,  and  the  difference 
equations  are  accurate  to  order  (AC)^,  (AS)^,  respectively.  The  areas  of 
the  cell  sides  are  the  same  as  before,  e.g.  , 


- X 


(48) 


'i+1 


INITIAL  CONDITIONS 


With  any  numerical  method,  the  closer  the  initial  guess  is  to  the 
answer,  the  faster  the  convergence.  In  Che  present  scheme  a simplified 
form  of  Slender  Body  Theory,  using  superimposed  solutions  from  the  linear- 
ized potential  equation,  is  employed  to  generate  initial  values  of  4)  and  p 
in  the  grid.  Each  centerline  segment  of  the  body  lying  In  a cell  is  assumed 
to  be  a constant  strength  source  filament.  The  strength  is  proportional 
to  Che  value  of  the  body  surface  alone  as  evaluated  at  the  midpoint  of  the 
segment.  The  initial  value  of  4>  at  any  point  is  then  the  sum  of  the  contri- 
butions from  each  segment. 

COMPARISON  WITH  EXPERIMENTAL  DATA 

Figures  7 through  10  present  comparisons  of  theory  with  experiment  for 
several  two-dimensional  and  three-dimensional  axlsymmetric  bodies.  The 
figures  have  experimental  and  theoretical  values  of  pressure  coefficient 
Cp  plotted  versus  non-dimensional  body  lengths. 

Figure  7 presents  data  for  a 6- percent- thick  parabolic  arc  airfoil  in 
two-dimensional  flow.  Figure  7(A)  has  a free  stream  mach  number  of  0.707, 
and  Figure  7(B)  has  M^  equal  to  0.817.  The  body  is  twenty  cells  long  and 
one  cell  high.  In  both  cases  the  flow  is  subcrltical,  and  the  critical 
pressure  coefficient  Cp*  is  listed  on  the  figure  and  shown  on  the  axis  where 
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O EXPERIMENT  (REFERENCE  16) 


(C).  Schematic  of  Body  Profile 


Figure  7.  Comparison  of  Numerical  Solution  Using  Purvis  Method 
and  Experiment  for  a 6- Percent  Parabolic  Arc  Airfoil 


2 


(Ci  Schematic  of  Double  Wedge  Airfoil 


Figure  8 


Comparison  of  Numerical  Solution  Using  Purvis  Method 
and  Experiment  for  Double  Wedge  Airfoil 


G EXPERIMENTAL  (REFERENCE  13) 
THEORETICAL  (PURVIS  METHOD) 


appropriate.  The  flow  is  close  to  reaching  critical  in  Figure  7(B).  The 
agreement  between  theory  and  experiment  is  excellent  in  both  cases. 

Figure  8 presents  the  same  type  of  data  for  two-dimensional  flow  over 
a double  wedge  with  a half-angle  of  4.3  degrees.  The  body  is  twenty  cells 
long  and  one  cell  high.  Note  that  even  with  this  coarse  grid  the  theory 
attempts  to  predict  the  Infinite  Cp  on  the  sharp  corner  at  the  center  of  the 
wedge.  In  Figure  8(A),  the  critical  pressure  coefficient  is  reached  near 
the  wedge  maximum  thickness.  The  agreement  between  theory  and  experiment 
deteriorates  rapidly  in  this  region.  The  reason  for  this  poor  agreement  is 
that  the  current  method  has  no  mechanism  such  as  upwind  differencing  to 
properly  account  for  the  zone  of  dependence  of  the  differential  equations 
and  to  enforce  the  directional  property  of  the  second  law  of  thermodynamics. 

Figure  9 presents  data  for  the  M117  standard  bomb  at  free  stream  mach 
numbers  of  0.5,  0.8,  and  0.9.  The  M117  is  a fairly  blunt  axlsymmetrlc  body 
consisting  of  a tangent  ogive  nose  section,  a cylindrical  midsection,  and  a 
straight-tapered  boattall.  The  body  length  is  nondlmensionalized  with 
respect  to  the  body  diameter  d.  In  Figure  9(A),  the  body  is  48  cells  long 
and  in  Figure  9(B)  the  body  is  40  cells  long;  in  both  cases  the  body  is 
5 cells  high.  The  effect  of  sting  location  on  the  pressure  distribution  is 
Illustrated  comparatively  in  the  two  figures.  The  theoretical  sting,  shown 
in  the  schematic  below  each  plot,  begins  at  different  locations  for  the  two 
cases,  4.8  body  diameters  aft  of  the  nose  in  Figure  9(A)  and  4.0  body 
diameters  aft  in  Figure  9(B).  The  predicted  and  measured  pressure  distribu- 
tions Indicate  that  the  sting  actually  begins  at  about  4.3  diameters.  The 
agreement  between  theory  and  experiment  is  excellent  in  both  cases,  even  at 
the  high  Mach  number,  and  particularly  over  the  nose  of  the  body.  In  Figure 
9(C),  the  critical  mach  number  has  been  exceeded  and  a shock  wave  has  formed 
at  an  axial  position  of  about  x/D  ■ 1.6.  As  expected,  the  shock  wave  and 
its  associated  pressure  jump  cannot  be  predicted  by  potential  theory. 

Figure  10  presents  data  for  three-dimensional  axlsymmetrlc  flow  over  a 
maximum  volume  bomb  at  a free  stream  mach  number  of  0.5.  The  bomb  consists 
of  a hemispherical  nose,  a straight  cylindrical  midsection,  and  a tangent 
ogive  faired  into  a straight  tapered  boattall.  This  configuration  is  a 
severe  test  of  the  theory  due  to  the  nose  bluntness  and  curvature.  Overall 
agreement  between  theory  and  experiment  is  excellent  for  this  case  particu- 
larly on  the  nose.  The  slight  differences  on  the  aft  section  of  the  body 
have  been  attributed  to  tunnel  flow  angularities.  This  particular  configur- 
ation best  Illustrates  the  versatility  and  accuracy  of  the  method. 

COMPUTATIONAL  PARAMETERS 

Further  comparisons  of  significant  computational  parameters  between  the 
various  cases  and  two  other  finite-difference  methods  are  given  in  Table  1. 
The  execution  times  for  the  bodies  in  Figures  7 through  10  varied  from 
3 to  10  minutes,  depending  on  the  body  shape  and  free  stream  mach  number. 

The  execution  times  generally  increased  as  the  free  stream  Mach  number  was 
increased.  The  longest  run  time  was  for  the  maximum  volume  bomb;  this  was 
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Figure  10.  Comparison  of  Nuiierlcal  Solution  Using  Purvis  Method 

and  Experlisent  for  16-Inch  Maximum  Volume  Hemispherical 
Nose  Bomb 


attributed  to  the  large  velocity  gradients  In  the  nose  region.  In  separate 
computer  runs,  the  execution  times  for  the  M117  bomb  decreased  from  ten  to 
twenty  percent  when  Slender  Body  Theory  was  used  for  Initial  conditions 
rather  than  starting  the  problem  from  free  stream  values. 

For  all  cases  shown,  the  grid  size  was  140  by  50  cells.  The  number  of 
cells  along  the  bodies  In  the  x-directlon  varied  from  20  to  68.  Separate 
computer  runs  showed  that  coordinate  stretching  was  unnecessary  as  long  as 
embedded  supersonic  regions  did  not  appear  In  the  flow. 

COMPARISON  WITH  OTHER  METHODS 

Computational  parameters  for  two  other  finite-difference  methods  are 
also  presented  at  the  bottom  of  Table  1.  Results  from  the  pioneering  work 
of  Murman-Cole  (Reference  16)  are  shown  for  the  same  6-percent  circular  arc 
airfoil  in  Figure  7.  The  execution  time  is  much  greater  even  though  their 
method  had  fewer  total  grid  points  (or  cells) . The  method  used  the 
transonic  small-perturbation  equation  rather  than  the  full  potential  equa- 
tion. Their  method  was  also  more  complex  and  used  nonlinear  coordinate 
stretching  in  both  directions,  circle-plane  mapping  and  special  treatment 
of  the  far-field  boundaries. 

Results  of  South  and  Jameson  (Reference  11)  are  given  for  a body 
similar  to  the  M117  bomb:  a tangent  ogive  nose,  cylindrical  midsection, 
and  straight  tapered  boattall.  No  specific  execution  time  for  this  case 
was  given  so  an  average  estimate  was  made  based  on  given  run  times  per 
grid  point  and  total  number  of  Iterations.  The  results  Include  successive 
mesh  halving,  although  they  state  chat  without  mesh  halving  the  run  times 
Increase  by  at  least  a factor  of  five  and  the  accuracy  decreases.  The 
method  used  three  different  coordinate  systems  on  the  various  sections  of 
Che  body.  It  was,  however,  designed  for  transonic  solutions  and  did  solve 
the  full  potential  equation. 

CONCLUSIONS 

The  present  method  is  computationally  fast,  accurate,  and  compares 
well  with  experiment.  Computation  times  are  better  than  other  methods  when 
all  factors  such  as  specific  computer,  number  of  grid  points,  and  special 
convergence  techniques  are  accounted  for.  The  computer  program  Is  simple 
and  Involves  no  complicated  numerical  schemes  or  empirical  techniques.  For 
all  cases  considered,  coordinate  stretching  was  found  to  be  unnecessary, 
although  the  method  Is  restricted  to  flows  with  no  embedded  supersonic 
regions.  No  problems  were  encountered  for  blunt  body  cases  or  near  r»0  for 
axlsymmetrlc  flows,  and  no  special  treatment  was  necessary  for  these  regions 
or  for  handling  the  boundary  conditions  on  body  surfaces. 

RECOMMENDATIONS 


The  present  method  has  proven  to  be  sufficiently  fast,  simple,  and 
accurate.  For  solving  large  problems  such  as  two-  and  three-body  inter- 
ference flows,  the  following  modifications  are  recommended: 
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(1)  Use  of  line  relaxation  to  reduce  storage  and  perhaps  accelerate 
solution  convergence. 

(2)  Incorporate  a successive  mesh  refinement  scheme,  such  as  that  of 
Reference  17  which  could  reduce  execution  time  by  a factor  of 
five. 

(3)  Use  a linearized  planar  wing-pylon  solution  for  wing-store  inter- 
ference. 

(4)  Use  previous  methods  such  as  Reference  6 for  initial  conditions 
and  far-field  boundary  conditions. 

(5)  Add  upwind  differencing  or  some  equivalent  scheme  to  allow 
extension  of  the  method  above  critical  mach  numbers. 


SECTION  IV 


0 


t 


DISCUSSION  OF  BOUNDARY  LAYER  STUDY 

INTRODUCTION 

As  can  be  seen  by  reading  the  literature,  no  one  to  date  has  solved 
the  three-dimensional  turbulent  boundary  layer  shock  Interaction  problem. 

The  basic  reason  stems  from  two  distinct  areas:  (1)  the  mathematical 
modeling  of  the  flow  field  and  (2)  the  viscosity  model  for  turbulent  flow. 

In  spite  of  these  problems,  however,  attempts  to  model  the  flow  In  the 
present  work  were  carried  out  using  a simple  eddy  viscosity  formulation. 

Most  of  the  work  was  done  In  generating  an  external  flow  field  and  coupling 
the  Invlscld  field  to  the  viscous  boundary  layer.  Because  of  computation 
time  and  storage  on  the  computer,  techniques  were  used  which  would 
require  minimal  time  and  storage.  For  this  reason,  solutions  of  the  Navler 
Stokes  equations  were  eliminated  from  consideration  at  the  outset.  It  was 
also  decided  that  (1)  the  general  Interference  potential  flow  solution 
already  developed  In  past  research  efforts  would  be  utilized  to  generate 
the  basic  Invlscld  solution;  (2)  three-dlmenslonallty  must  be  considered; 
and  (3)  compressible  turbulent  flow  must  be  Incorporated  Into  the  solution. 
Because  of  these  restrictions,  the  boundary  layer  equations  approach  was 
used  In  the  formulation  of  the  viscous  flow  region.  In  this  region,  the 
basic  equations  were  taken  from  the  work  of  Sasman  and  Crescl  (Reference  18) . 

MATHEMATICAL  APPROACH 

The  Idea  In  the  present  work  was  to  develop  a technique  which  would 
couple  the  present  three-dimensional  potential  flow  solution  of  the  three 
stores  mounted  on  a triple  ejector  rack  to  a ttro-dlmenslonal  boundary 
layer  formulation.  This  concept  Is  feasible  If  discrete  streamlines 
located  at  the  edge  of  the  boundary  layer  could  be  Identified  as  shown  In 
Figure  11.  Since  the  boundary  layer  Is  generally  thin,  the  layer  of 
fluid  Immediately  below  the  streamline,  l.e.,  between  the  streamline  and 
the  actual  body  surface,  can  be  considered  to  be  two-dimensional  in  nature 
In  which  the  two-dimensional  boundary  layer  equations  apply.  That  la, 
even  though  the  streamline  may  curve  around  the  body  and  the  plane  defined 
by  the  streamline  and  the  shadow  of  the  streamline  on  the  body  Is  a curved 
surface,  the  thickness  Is  much  less  than  the  body  characteristic  dimension 
and  thus  may  be  approximated  by  two-dimensional  boundary  layer  equations. 
Hence,  a clear  accurate  definition  of  a streamline  near  the  body  surface 
Is  essential  to  the  problem  solution. 

The  mechanics  of  generating  this  streamline  is  rather  simple.  First, 
the  potential  flow  solution  Is  generated  for  the  system  of  bodies  chosen 
to  be  analyzed.  Then  a point  In  apace  (xo,yo,ZQ)  Is  chosen  near  the  body 
surface  and  near  the  nose  of  the  body.  AC  this  point,  using  the  potential 
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solution  already  determined,  the  velocity  components  u' , v',  and  w'  are 
determined.  Adding  the  free  stream  components  appropriately,  we  then  have 


u » cos  ot  + u'  (49) 

V - sin  a + v'  (50) 


w 


w' 


(51) 


where  a Is  the  body  angle  of  attack.  Then  It  Is  assumed  that  over  a very 
small  but  finite  time  Interval  the  new  positional  coordinates  of  the 
streamline  may  be  approximated  by 


t 


} 


and 
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X + uAt , 
o 

y + vAt , 
o 

z + wAt, 
o 


(52) 


Computation  of  the  new  fluid  velocity  at  the  new  location  Is  then 
carried  out,  and  the  process  Is  repeated.  In  this  manner  a matrix  of 
streamline  coordinate  locations  Is  defined.  This  method  suffers  from  one 
major  drawback.  Any  error  Induced  in  the  coordinate  location  at  any  up- 
stream point  Is  carried  to  the  next  location  and  accumulates.  However, 
large  error  buildup  can  be  avoided  by  the  use  of  small  time  intervals.  At. 
The  coordinates  of  these  points  on  the  streamline  are  also  used  to  determine 
the  step  Interval  for  the  boundary  layer  computations,  and  the  velocity  at 
the  edge  of  the  boundary  layer  Is  always  known  as  well  as  the  external 
free  stream  pressure. 

This  section  of  the  program  and  analysis  proved  to  be  rather  success- 
ful In  that  discrete  streamlines  beginning  anywhere  near  the  body  surface 
were  Identified.  Problems  were  encountered  near  the  aft  end  of  the  body 
where  the  streamline  deviated  from  the  body  and  trailed  off  In  the  free- 
stream  direction.  This  problem  was  attributed  to  the  characteristic  nature 
of  point  source  solutions  for  the  potential  flow  field  in  which  the  induced 
velocity  component  u tends  to  oscillate. 

Boundary  layer  C'jmpu tat ions  using  this  scheme  were  carried  out  with 
some  degree  of  success.  However,  the  interaction  problem  using  this 
approach  has  not  been  solved. 
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CONCLUSIONS 


After  several  computer  runs,  more  literature  studies,  and  much  thought. 
It  Is  believed  that  the  approach  taken  here  will  not  yield  the  desired 
results.  This  conclusion  Is  based  on  the  following  observations: 

1.  The  shock  boundary  layer  Interaction  Is  strongly  dependent  on  Its 
three-dimensional  geometry. 

2.  In  the  boundary  layer  equations  approach,  the  boundary  layer  only 
reacts  to  Its  external  flow  field;  hence  this  field  must  be  very  accurate  If 
accurate  results  are  obtained.  Provisions  using  this  scheme  must  be  made 
for  transmission  of  disturbances  upstream  In  the  boundary  layer  or  In 
regions  where  M < 1.0  (elliptic  equations  In  the  potential  field)  and  must 
obey  the  hyperbolic  equations  for  regions  where  M > 1.0. 

3.  At  present,  the  Navler  Stokes  solution  approach  requires  too  much 
computer  time  and  storage. 

4.  Turbulent  mixing  models  are  Indequate  to  model  the  shear  stresses 
in  the  interaction  region. 

5.  It  appears  that  much  more  work  In  viscous  aerodynamics  will  be 
required  before  adequate  results  will  be  obtained  by  pure  analytical 
methods . 
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SECTION  V 


DISCUSSION  OF  FINITE  ELEMENTS  STUDY 


INTRODUCTION 


The  analysis  of  subsonic  flow  requires  the  solution  of  partial 
differential  equations  of  the  elliptic  type.  The  solution  of  these 
equations  can  be  obtained  by  the  traditional  Finite  Difference  Methods, 
as  mentioned  in  Section  111.  One  of  the  Inherent  difficulties  with  the 
Finite  Difference  Method,  which  is  not  significant  in  the  finite  element 
method,  is  that  special  differencing  techniques  must  be  incorporated  to 
handle  Irregular  or  curved  boundaries.  Since  the  difference  formulas  used 
for  these  boundaries  are  usually  of  lower-order  accuracy  than  the  formulas 
used  for  the  rest  of  the  flow  field  and  since  in  subsonic  flow  disturbances 
at  any  point  in  the  flow  region  are  felt  throughout  the  flow  region,  the 
overall  accuracy  of  the  solution  is  reduced.  This  inexact  treatment  of  the 
boundaries  is  a limitation  to  the  Finite  Difference  Method  although  very 
accurate  solutions  may  be  obtained  by  using  a large  number  of  grid  points 
or  by  using  higher  order  difference  formulas.  However,  the  incorporation 
of  more  accurate  formulas  in  the  Finite  Difference  Method  usually  compli- 
cates the  treatment  of  the  boundaries. 

In  the  last  decade  much  effort  has  been  directed  toward  the  application 
of  the  Finite  Element  Method  to  fluid  flow  problems.  Although  it  was 
originally  developed  in  the  aircraft  industry  for  solution  of  complex 
structural  problems,  the  Finite  Element  Method  is  readily  applicable  to 
all  types  of  boundary  value  problems  (see  Reference  19) . An  application  of 
this  method  to  Incompressible  Invlscld  potential  flows  is  given  by  Martin 
(Reference  19) . Zienkiewlcz  (Reference  15)  gives  an  application  to 
Incompt esslble  viscous  flow  problems.  Huebner  (Reference  20)  devotes  an 
entire  chapter  to  various  methods  which  have  been  successfully  used  to 
solve  both  viscous  incompressible  and  Invlscld  compressible  flows.  Ref- 
erence 21  presents  numerous  analyses  of  fluid  flow  problems  which  have  been 
solved  by  the  Finite  Element  Method. 


^ A particular  formulation  given  by  Heubner  (Reference  20)  has  been  used 

' by  a number  of  researchers  to  obtain  solutions  to  Invlscld  compressible 

t flows.  This  formulation,  in  a slightly  modified  form,  has  been  chosen  in 

i ■ this  study  to  obtain  solutions  to  two-dimensional  flow  over  a circular  arc 

I airfoil  for  mach  numbers  where  compressibility  effects  are 

> governing  equation  and  the  basic  iteration  are  the  same  as 

1 . Section  III  of  this  report.  The  formulation  and  numerical 

I ever,  are  significantly  different. 


significant.  The 
those  used  in 
technique,  how- 


GOVERNING  EQUATIONS 


The  governing  equations  in  this  study  are  the  same  as  used  in 
Section  III.  Repeated  here  for  continuity,  the  conservation  equations 
for  steady,  adiabatic,  inviscid  flow  of  an  ideal  gas  may  be  combined  to 
give  the  equations 


V*(pV<(.)  - 0 (55) 

and  ^ 

p - [1  + (^)  (1  - (56) 

where  (|i  is  the  nondlmenslonal  velocity  potential,  and  p is  the  nondimen- 
sional  density. 

SOLUTION  PROCEDURE 

Equations  (55)  and  (56)  are  the  basis  of  a rapidly  convergent 
iterative  method  for  solution  of  compressible  flow  problems.  The  first 
step  involves  assuming  p • 1 throughout  the  flow  field  as  an  Initial 
estimate.  Then  Equation  (55)  is  used  to  derive  a system  of  linear  alge- 
braic equations  which  may  be  solved  for  values  of  velocity  potential  (Ji. 

The  resulting  values  of  <p  are  then  numerically  differentiated  to  obtain 
values  of  the  velocity  components  which  are  used  in  Equation  (56)  to 
compute  new  values  of  density  which,  in  turn,  are  used  in  Equation  (56)  and 
the  process  is  repeated.  The  iteration  is  continued  until  the  densities  no 
longer  change  from  one  iteration  to  the  next.  The  pressures  are  then 
computed  using  Equation  (57) . 

C - — {[1+^m2  (1  - V4.-V^)]  -1}  (57) 

OP 

which  is  the  usual  expression  for  pressure  coefficient  for  compressible  flow. 
FINITE  ELEMENT  FORMULATION 

The  flow  field  in  which  the  solution  is  obtained  is  an  infinite  domain 
surrounding  a two-dimensional  6-percent  circular  arc  airfoil.  The  far-field 
domain  is  approximated  by  placing  the  computational  boundaries  three  body 
lengths  away  from  the  body  in  all  directions  and  forms  a rectangular  region 
enclosing  the  body.  The  far-fleld  boundary  condition  used  Is  that  of 
uniform  flow  in  the  x direction  and  expressed  in  terms  of  the  velocity 
potential  is  given  by  Equation  (58) 


The  body  boundary  condition  Is  the  usual  flow  tangency  condition  expressed 
by  Equation  (59) 


cV4i*n  - 0 (59) 

where  n Is  a unit  vector  normal  to  the  body  surface. 

‘ For  computational  purposes,  these  mixed  boundary  conditions  are 

formulated  as  given  by  Equations  (60)  and  (61) . 

pV<|i*n  - f(x,y,z)  on  Sj  (60) 

and 

't'  “ g(x,y,z)  on  (61) 

where  and  $2  are  portions  of  the  boundary  and  f and  g are  prescribed 
functions  of  position  (see  Figure  12) . 

The  finite  element  solution  technique  Is  based  on  a variational 
principle  derived  from  the  calculus  of  variations  which,  when  expressed  In 
terms  of  the  velocity  potential,  minimizes  the  functional 

K-li)  - 7 J pV(t.-7(j.  dO  - / (}.f  dS,  (62) 

fi  ®2 

and  satisfies  Equations  (55),  (56),  and  (61). 

' To  apply  this  method,  the  rectangular  region  Is  subdivided  Into 

triangles  and  the  values  of  <p  are  Interpolated  over  the  area  of  the 
r triangle  by  means  of  a set  of  Lagranglan  Interpolation  functions.  These 

Interpolation  formulas  are  most  conveniently  represented  In  terms  of  the 
I natural  coordinates  (sometimes  called  area  coordinates)  defined  by  the 

equations 


X - L,  x,  + L,  X-  + L-  x (63a) 


where  (x  ,y  , (X3,y3)  are  the  vertices  of  the  triangle  In  plain 

rectangular'^^cooralnates.  Solving  Equations  (63a,  b,  c)  for  L^,  L.,  and 
yields  the  general  form 

Lf  - Hi  + bf  X + Cf  y,  i - 1,  2,  3 (64) 

where  a^,  b^,  and  are  constants. 

For  linear  Interpolation  between  vertices  of  the  triangle,  the  Inter- 
polation functions  are  given  by 


N 


1 


(65) 


For  second-order  Interpolation  they  are  given  by 


N 

N 

N 

N 

N 

N 


1 

2 

3 

4 

5 

6 


Ll  - L2(L3+L^) 


l2  - L3(L,+L2) 


4L3L2 


4L2L3 


4L3L1 


(66) 


where  the  node  numbers  1 to  6 are  shown  In  Figure  13. 

The  nodal  values  of  are  used  to  define  the  variation  of  ^ over  a 
triangular  element  by  the  following  expression 


?"i*i 


(67) 


After  the  Interpolation  functions  are  chosen,  they  are  substituted  into 
Equation  (62)  and  then  che  mlnisdzatlon  process  Is  carried  out  on  an 
element-by-element  basis.  This  is  done  by  requiring  that 
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34> . 


■ 0 5 j “ 1»  2, 


(68) 


where  r ■>  3 for  linear  Interpolation  and  r « 6 for  second-order  Interpola- 
tion. 


After  performing  the  indicated  differentiations,  there  is  obtained 
for  each  element  a set  of  equations  of  the  form 


[K]{ij)}  - {R} 


(69) 


where  is  an  r x r matrix  of  Influence  coefficients  relating  the  nodal 

variables  4)^ , which  is  an  r x 1 column  vector  of  nodal  values  of  velocity 

potential,  and  is  an  r x 1 column  vector  of  forcing  functions.  For  this 

case  ” 0 since  f(x,y,z)  ■ 0 on  the  body  surface.  The  k-matrix  is 

obtained  by  first  deriving  the  element  matrix  equations  and  then  assembling 
these  elements  into  a system  matrix  containing  the  summation  of  all  of  the 
element  contributions  to  each  nodal  variable. 


The  resulting  system  of  equations  must  be  modified  to  incorporate  the 
nodal  boundary  conditions  as  defined  in  Equation  (61).  To  apply  a known 
value  of  4>  at  the  1th  node,  the  1th  system  equation  is  cleared  to  zero. 

Then  the  diagonal  entry  is  set  to  one  and  the  known  value  placed  on  the 
right-hand  side  of  the  system  equations  in  the  ith  position.  This  forces 
the  solution  of  the  system  to  obey  the  specified  boundary  conditions  at  the 
ith  node. 

After  application  of  the  boundary  conditions  to  the  system,  the  equa- 
tions are  solved  by  means  of  a banded  Gauss  elimination  procedure.  The 
band  width  of  the  system  matrix  is  highly  dependent  upon  the  nodal  nundier- 
Ing  scheme  and  is  dependent  solely  upon  the  maximum  difference  in  nodal 
numbers  within  any  given  element  taken  over  all  the  elements.  It  was  found 
that  the  run  time  and  core  storage  requirements  were  significantly  reduced 
with  no  loss  of  accuracy  by  taking  into  account  the  banded  nature  of  the 
system  matrix.  All  of  the  calculations  were  performed  on  an  IBM  370/158 
computer  in  single  precision  mode. 

After  solving  for  new  values  of  the  results  were  numerically 
differentiated  by  means  of  the  Interpolation  functions.  For  exaiq>le. 


then  the  resulting  velocities  are  used  in  Equation  (56)  to  obtain  new  values 
of  density  at  every  point.  The  process  was  then  repeated  until  convergence 
was  attained. 
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After  the  velocities  were  obtained  from  the  final  Iteration,  the 
pressure  coefficients  were  obtained  from  Equation  (57) 

RESULTS  AND  CONCLUSIONS 

For  the  case  of  the  6-percent-thlck  circular  arc  airfoil  at  M “ 0.707, 
the  results  are  shown  in  Figure  14.  Figure  14(A)  shows  the  results  for  the 
second-order  interpolation  functions.  The  fiTst-order  solution  appears  to 
be  very  smooth  and  follows  the  experimental  trends  well  even  though  the 
finite  element  grid  is  quite  coarse.  It  would  be  expected  that  this  solu- 
tion would  approach  the  true  solution  asymptotically  as  the  grid  is  refined. 
This,  of  course,  would  require  a much  greater  computer  core  and  time. 

An  alternate  method  for  Improving  the  accuracy  of  the  solution  is  to 
use  second-order  interpolation  functions.  This  technique  was  implemented 
with  the  results  shown  on  Figure  14(6). 
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SECTION  VI 


SUMMARY  OF  RESULTS  WITH  COMMENTS  AND  RECOMMENDATIONS 


SUMMARY  OF  RESULTS  WITH  COMMENTS 

The  Raylelgh-Janzen  solution  was  generated  through  the  second-order 
terms  for  small  disturbances  and  through  the  first  order  for  a modified 
small  disturbance  formulation.  Although  It  had  been  anticipated  that  this 
type  of  solution  would  be  more  applicable  to  flow  conditions  near  the 
critical  mach  number,  the  resulting  pressure  predictions  are  no  better 
than  those  obtained  by  using  Gothert's  Rule.  In  addition.  It  was  found 
that  although  small  disturbances  are  assumed,  both  methods  predict  quite 
well  the  pressures  In  the  vicinity  of  stagnation  points.  The  computer  time 
and  core  requirements  are  about  the  same  for  the  Raylelgh-Janzen  solutions 
as  for  the  solutions  using  Gothert's  Rule,  approximately  three  to  five 
seconds  on  a CDC  6600  computer. 

Preliminary  work  Indicates  that  the  Iterative,  numerical,  successive 
approximation  method  outlined  In  Section  II  (based  on  the  Raylelgh-Janzen 
method)  may  give  a very  fast  way  to  correct  any  Incompressible  solution 
(zeroth  solution)  for  compressibility  effects  up  to  the  critical  mach 
nuiid>er. 

The  Purvis  cell  method  is  a numerical  technique  based  on  an  Integral 
form  of  the  conservation  equations  for  Irrotatlonal  flow.  Even  though  this 
solution  Is  Iterative  In  nature,  the  computational  time  Is  an  order  of 
magnitude  faster  than  finite  difference  methods.  An  additional  and 
significant  reduction  In  time  can  be  obtained  by  using  Slender  Body  Theory 
for  predicting  the  initial  estimate;  however,  this  method  still  takes  an 
order  of  magnitude  longer  than  Slender  Body  Theory.  The  agreement  with 
experimental  results  Is  excellent  for  mach  numbers  up  to  and  slightly 
greater  than  the  critical  value.  This  method  Is  applicable  to  any  general 
body  shape  and,  therefore,  should  be  quite  useful  for  future  applications. 

The  state-of-the-art  of  viscous  flow  analysis  does  not  seem  to  be 
suitable  for  use  on  complicated  flow  fields  such  as  associated  with  triple 
ejector  racks,  at  least  within  the  constraints  of  limited  computer  core  and 
computational  time. 

The  finite  element  method  seems  to  work  quite  well  for  potential  and 
some  viscous.  Incompressible  flow  cases.  The  application  required  for 
mutual  Interference  - compressible,  viscous  flow,  however,  has  not  been 
successfully  formulated. 
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RECOMMENDATIONS 


Several  approaches  Investigated  In  this  work  warrant  additional  study 
as  suggested  below; 

1.  The  successive  numerical  approximation  outlined  in  Section  II 
should  be  Investigated.  This  technique  could  be  used  to  correct  any 
Incompressible  flow  solution  for  compressibility  effects. 

2.  The  Purvis  cell  method  should  be  expanded  to  Include  the  general 
three-dimensional  problem.  Further,  It  Is  recommended  that  a coordinate 
transformation  to  rectangular  arrays  be  considered  In  order  to  simplify 
the  utilization  of  this  method  for  arbitrary  curved  bodies. 

3.  The  potential  flow  solution  should  be  modified  so  as  to  Include 
shock  jump  effects. 

It  should  be  noted  that,  although  the  Purvis  method  takes  an  order  of 
magnitude  longer  computer  time  than  the  Raylelgh-Janzen  or  Gothert  solutions 
the  velocities  on  the  body  surface  are  computed  as  part  of  the  solution. 
Thus,  the  complete  pressure  distribution  Is  available  as  part  of  the  solu- 
tion. On  the  other  hand,  the  times  quoted  for  the  Raylelgh-Janzen  and 
Gothert  solutions  are  for  Che  basic  solution  giving  the  strength  distribu- 
tion of  Che  source  singularities.  VAien  these  basic  solutions  are  used  to 
compute  the  complete  body  pressure  distribution,  the  computational  times  are 
the  same  order  of  magnitude  as  for  the  Purvis  method. 
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APPENDIX  A 

RAYLEIGH- JANZEN  SOLUTION  FOR  SLENDER  BODIES 

The  velocity  potential  equation  can  be  written  in  a form  useful  in 
this  study  to  yield 


where 


a2v2^  - - (:^)  V^.Va 


a2  - a2tl+(^)M2(l-7<|,.7$)] 


(A-1) 


(A-2) 


Expanding  Equation  (A-1)  in  terms  of  the  nondlmenslonal,  perturbation 
velocity  potential  yields 

[1-(^)m2(2^^-H2-N>2-H>2)]v2^  . 


f i III  **l**l* 


(A-3) 


Using  e ■ m2  as  the  perturbation  parameter,  an  assumed  series  solution 
of  the  form 


.(0) 


(>(x,y,z;e)  - ♦""(x,y,z)  + (x.y.z)  + 


+ e2^^2)(x,y,z)  + ... 


(A-4) 


will  be  considered  as  follows. 


The  assumed  solution  can  be  substituted  into  Equation  (A-3)  and,  in 
principle,  each  successive  order  of  magnitude  term  can  be  solved.  For 
small  perturbations,  however.  Equation  (A-3)  can  be  simplified  by  neglecting 
velocity  squared  terms  In  comparison  to  first-order  terms  and  Equation 
(A-3)  reduces  to 

II  - (Y-l)  >t  ♦,!  V** . ,ei<v»  I;  Ill  ♦.  (A-5) 


4S 


The  right-hand  aide  of  Equation  (A-5)  can  now  be  rewritten  Into  the 
equivalent  form  given  by  Equation  (A-6) 

4i2  + (j,2  + ^2 

[1  - (y-1)  m2  v2^  - m2  1^  ] (A-6) 

Again  neglecting  second-order  velocity  terms  in  comparison  to  the 
first-order  terms  yields 

[1  - (Y-1)  m2  V2^  . m2  (A-7) 

The  body  boundary  conditions  for  axlsymmetrlc  bodies  at  zero  angle  of 
attack  can  be  expressed  as 


dr 


cos  6 + ® ^ 


(A-8) 


body 


Substituting  Equation  (A-4)  into  Equations  (A-7)  and  (A-8),  and  equating 
terms  in  powers  of  e gives  the  first  three  orders  of  magnitude  equations  as 
follows: 


ZEROTH  ORDER 


72  (ji 


(0) 


(0) 


FIRST  ORDER 


cos 


e + 511.  6 - (!+♦*“*)  g 


. ♦<»> 

^ ’^XX 


body 


cos  e + sin  0 - ^ 


body 


SECOND  ORDER 


,2,<«  . ♦«»  *«»  + 


COS 

y 


0 + > sin  0 - 


body 


(A-9a) 

(A-9b) 

(A-lOa) 

(A-lOb) 

(A-lla) 

(A-llb) 
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SOLUTION  FOR  AXISYMMETRIC  BODIES  (a-0) 


i 

i 


For  axisynmetrlc  bodies  at  zero  angle  of  attack,  the  solution.  In- 
cluding the  image  system  for  mutual  Interference,  can  be  expressed  as  the 
sum  of  the  source  singularities  as  given  by  Equation  (A-12)  for  the  zeroth 
solution. 


-C 


.(0) 


1 ’*1 


where 


-Y(x-x^)2  + (y-y^)2  + (z-z^)2 


(A-12) 


(A-13) 


and  the  source  strengths,  C , are  determined  so  as  to  satisfy  the  body 
boundary  conditions.  ”l 

The  first-order  solution  which  satisfies  Equation  (A-10)  can  be 
expressed  as 

(1)  r 1 V 

^ - I Hr  2 ^ 

i i i R^^ 

%diere  the  source  strengths,  C.  , are  determined  by  satisfying  the  flrst- 
orUer  boundary  conditions.  j 

After  using  Equations  (A-12)  and  (A-I4)  In  (A-11),  the  second-order 
solution  can  be  found  to  be 

C 


-C, 


♦ 


(2) 


I 

k 


A.iA  L_ 

c Rl 


C,  (x-xj^  (x-x^)‘ 


^ J b] 


+ (Y-1)  E — 


C (x-x  ) 
o,  1 


“I  »i 

C (x-x.) 2 

y l_i 

^2  1 

1 rJ 


(A-15) 


The  source  strengths,  C2  • must  be  found  so  as  to  satisfy  the  second- 
order  boundary  conditions.  ^ 
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The  constants,  C^j  and  C2|^,  are  the  source  singularity  strengths  for 

the  complimentary  solutions  of  Equations  (A-10)  and  (A-11),  respectively. 
It  Is  possible,  however,  that  In  addition  to  (or  rather  than)  the  sources, 
doublet  singularities  niay  be  required  in  order  to  satisfy  the  boundary 
conditions.  If  that  is  the  case.  Equations  (A- 14)  and  (A-15)  will  require 
terms  such  as  given  by  Equation  (A-16) . 


doublets 


+ I _n D_ 


(A-16) 


m 


m 


The  three  terms  In  Equation  (A-16)  are  for  doublets  with  their  axis 
along  the  x-,  y-,  and  z-dlrections,  respectively. 

AXISYMMETRIC  BODY  AT  ANGLE  OF  ATTACK,  a 

For  axlsymmetrlc  bodies  at  angle  of  attack,  the  solution  can  be  ex- 
pressed as  the  sum  of  source  singularities  plus  the  sum  of  doublet  singular- 
ities which  are  needed  in  order  to  satisfy  the  body  boundary  conditions  due 
to  the  cross-flow.  In  this  case,  the  body  boundary  conditions  can  be 
expressed  by  Equation  (A-17). 


dr 

A COS  0 + (iji  + sin  a)  sin  0 - ((J>  + cos  o)  — 

V z X dx 


body 


(A-17) 


The  order  of  magnitude  boundary  conditions  obtained  by  substituting 
Equation  (A-4)  Into  Equation  (A-17)  becomes,  for  the  first  three  terms. 


.(♦W) 

+ sin  a) 

sin 

0 - 

(,(0) 

+ 

(A-18) 

z 

X 

body 

,<l>cos 

y 

e . 

sin 

6 “ 

X 

dr  1 
dx 

body 

(A-19) 

cos 

e * 

sin 

0 - 

♦x"' 

dr 

(A-20) 

y 

ex 

body 

The  zeroth  solution  for  this  case  can  be  expressed  as  Equation  ^,-.l) 

-C 


u (z-z.) 

,(0)  . j ^ . J -X_J- 
11  J 


'A-21) 


where  the  singularity  strengths  must  be  found  so  as  to  satisfy  Equation 
(A-18).  Note,  however,  that  the  ?eroth  order  solutions  can  be  obtained  by 
solving  the  axisymsietrlc-f low  problem  and  the  cross-flow  problem.  Splitting 
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i 


the  body  boundary  conditions  Into  the  axial  and  cross-flow  components 
yields  for  the  axial  flow  problem. 


cos  0 + 


sin  0 


(f 


(0) 


+ cos  o) 


dr 

dx 


body 


(A-22) 


where  only  the  source  singularities  are  considered. 

For  the  cross-flow,  the  body  boundary  conditions  are  expressed  by 
Equation  (A-23)  where  only  the  doublet  singularities  are  considered. 


cos  0 +((j> 


(0) 


+ sin  a) 


sin  0 


^x  dx 


body 


(A-23) 


The  axial-flow  solution  remains  the  same  as  that  found  by  Equation 
(A-11)  except  that  the  boundary  condition  equation  must  be  corrected  by 
using  the  axial  component  of  the  free-stream  speed  as  given  by  Equation 
(A-22) . 


The  cross-flow  solution  is  now  given  by  Equation  (A-24)  where  Equation 
(A-23)  is  used  to  determine  the  doublet  singularities. 


(0) 


M (z-2  ) 

O,  1 


cross-flow 


(A-24) 


The  first-order  equations  can  also  be  solved  by  considering  the  axial- 
flow  and  the  cross-flow  solutions.  As  in  the  zeroth  solution,  the  axial- 
flow  solution  is  the  same  form  as  given  by  Equation  (A-14). 


The  governing  equations  for  the  cross-flow  problem  are  found  by 
substituting  Equation  (A-24)  into  Equation  (A-lOa)  and  yields 


cross-flow 


I 

k 


“k 


r5 


(A-26) 


The  second-order  equation  [see  Equation  (A-lla) ] contains  nonlinear 
terms  and,  therefore,  cannot  be  solved  by  separating  the  axial-flow  and 
cross-flow  solutions.  Hence,  the  governing  equations  for  this  case  is 


I 

» 

I . 

! -■ 
I 
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The  general  solution  for  the  second-order  term  can  be  expressed  as 


-C  n,  (z-z.) 

1 Rj  1 R] 


C (x-xj  3p^  (z-z  )(x-x  ) 
o.  i o.  J J 

(r-1)  ( I ♦ I — ; ) » 


i Ri 


»SC  (x-x.)2  u (x-x  )2(z-z  ) 

(V  ‘’i  — L 

I R?  ' j 


C (x-x  )2 

+ 1 y 

^1  R? 


C (x-x,)2 

*1  II  -^1 — 

^ i rJ 


fl^l ' 

^ i rJ 


+ I I 

k 


(z-Zj^)(x-x^^) 


■ 2 *^1, 


k Rj 


p (z-z.)(x-x  )•*  u (z-z  )(X-X  )- 

25'  »;  * «; 


(A-28) 


APPENDIX  B 


MODIFIED  SLENDER  BODY  SOLUTION 


The  governing  partial  differential  equation  which  properly  accounts 
for  stagnation  conditions,  small  velocity  perturbatlona  away  from 
stagnation  points,  and  conditions  up  to  the  critical  mach  number  [given 
in  Section  11  as  Equation  (20)]  is  given  here  as  Equation  (B-1) . 

(|>^  + (B-1) 

Equation  (B-1)  in  a slightly  more  convenient  form  becomes 

= M^[l  +(y-1)  (B-2) 


The  solution  of  Equation  (B-2)  can  now  be  generated  by  the  Rayleigh' 
Janzen  method  as  follows : 


Assuming  the  solution 


yields 

<»)  ; . 0 

(2)  . ,2,(2)  . * (»  + 7*«»-V*<2>  + V,<2).V,‘») 

* ^ XX  X X 

+ (Y-1) 


(B-3) 


(B-4a) 

(B-4b) 

(B-4c) 


S2 


A 


\ 

\ 


The  body  boundary  conditions  are: 


(0)  ; cos  6 + + sin  o)  sin  6 


(♦<“>  + cos  (B-5S) 


(1)  ; 


(2)  ; 


cos  e + sin  e 
y z 

♦ cos  0 + sin  e 
y * 


^(1)  (^) 

'dx^body 

^(2) 

^x  'dx^body 


(B-5b) 


(B5c) 


The  zeroth  solution  expressed  in  the  usual  manner  as  the  sum  of  sources 
and  doublet  yields 


♦«>  - 1 (♦<»>)  + I (♦<»>) 

I 'i  j j 

For  the  first-order  solution,  we  need  to  satisfy  the  equation 

y T 

Consider  first  the  following  expression 


(,(0),(0)) 


(B-6) 


(B-7) 


(B-8) 


which  can  be  expanded  to  give 


,2(, (»),(»))  . ,<«)  ,2,«»  + ,(0)  ,2,(0)  , 2„(0)  .,,,(0)  (M) 


. . 0 


Therefore 


,2(,(0),(0))  . 2,,'0'.,,« 


.(0).„JO) 


For  the  - term,  it  can  be  shown  that 

• mrme  * 


(B-10) 


(B-11) 


L 


(B-12) 


2 ^ 9 

1 rJ 


'^S  2 ^ 9 

«1  J r2 


Thus,  the  particular  solution  for  sources  and  doublets  can  be  written  as 
(x-x.)^  -1  _ fn\  ^ f fn'k  fO)  , 


!♦<») 


. i J ♦«»  . I I ♦<<»  ^ i I /«  (B-13) 


where  the  zeroth  solution  is 


♦"»- 1 * I C’ 

i ®i  j J 


♦<«>  - 1 /'’>  * I 

- 1 1 ■s.j 

The  complementary  solution  can  be  written  as 


„(»,  . I ♦<»  ^ J ,(1) 

1 i j j 

and  the  first-order  solution  becomes 


(B-U) 


(B-15) 


(B-16) 


(B-17) 
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